* Paddy Carter, Nicolas Van de Sijpe, and Raphael Calel, 13-June-2020
* Define program elusive used to run Monte Carlo simulations


/* This do file is called up by MCcountry.do so there is no preamble here. 
We use gtools to speed up the MCs run by elusive: https://github.com/mcaceresb/stata-gtools */


*************************************************************************
*************************************************************************
* ELUSIVE PROGRAM
*************************************************************************
*************************************************************************
/* Program key: */
/* (nC number of countries) (t years) (nD number of DFIs) (phi DFI country preferences weight) (drift in DFI budget) 
(sigma st. dev. in dfi budget) (dfbudg t=1 DFI budget as share of total opportunities)
(tvweight weight on time-varying part of DFI preferences: always zero in results reported in the paper but kept it in the code) 
(mu mean project characteristics of three country types) (musd st. dev. for project characteristics of three country types) 
(ersd st. dev. of expected returns noise added to pc) (dfilo lower bound eligible dfi set) (dfihi upper bound dfi set) 
(psmin min. exp return for private sector) (nI number of investment opportunities per country-year) (noise error added to pc after collapsing to country-year) 
estimator(estimator options are OLSFE, GMM, IV) (mechanism DFI investment selection mechanism: options are random, pclow)
transition(transition of country types over time; options are default, fixedtypes, oneperiodchange, fewertransitions) 
btot(to get billions to trillions effect: add (btot-1)*DFIinv to private investment) 
extraGMM(yes if want to run extra GMM models, yes only in core experiments, no for all the other experiments)
dfibudget(default or capped: capped constrains global DFI budget) */


capture program drop elusive
program define elusive, rclass
	version 15.1
	syntax [, nC(integer 90) t(integer 20) nD(integer 1) phi(real 0) drift(real 0) sigma(real 0.15) /// 
	dfbudg(real 0.1) tvweight(real 0) mu1(real 0) mu2(real 2) mu3(real 4) /// 
	musd1(real 1) musd2(real 1) musd3(real 1) ersd(real 1) dfilo(real 2) dfihi(real 1000) psmin(real 2) /// 
	nI(real 50) noise(real 0.5) estimator(string) mechanism(string) transition(string) btot(real 1) extraGMM(string) dfibudget(string)]

	
drop _all


* private sector budget (always exceeds available investment opportunities): 
local pb = `nI'*`nC'
* initial (year 1) value for the total DFI budget; defined as a fraction (dfbudg) of total number of annual investments:
local db = `dfbudg'*(`nI'*`nC')


/* Set up data: */	
* country panel:
local size = `nC'*`t'
set obs `size'
* number the countries:
gen country = ceil(_n/`t') 
* periods:
sort country	
by country: generate year = _n
xtset country year


**************************************************************************
* Generate weights to divide investments between multiple DFIs
**************************************************************************
* generate individual DFI country preferences that sum to 1 
forvalues j = 1/`nD' {
	sort year country
	gen double dfipref`j' = runiform() in 1/`nC' // generate DFI preferences for countries in first year only
	gegen double temp`j' = max(dfipref`j'), by(country) // this fills in preferences for all years
	drop dfipref`j'
	gen double tvtemp`j' = runiform() // time-varying part of DFI-country preferences. In the results we report in the paper, tvweight is always zero so this part drops out.
	gen double dfipref`j' = (1-`tvweight')*temp`j' + `tvweight'*tvtemp`j'
	drop temp`j' tvtemp`j'
}
ds dfipref*
egen double sum = rowtotal(`r(varlist)')
forvalues j = 1/`nD' {
	replace dfipref`j' = dfipref`j'/sum
}
drop sum


* generate individual DFI shares of each DFI investment in each country 
* start with weighted preferences: weight = pref^phi
* turn into shares using share_i = weight_i/(sum weights) 
forvalues j = 1/`nD' {
	gen double dfiweight`j' = dfipref`j'^`phi'
}
ds dfiweight*
egen double sum = rowtotal(`r(varlist)')
forvalues j = 1/`nD' {
	gen double dfishare`j' = dfiweight`j'/sum
}
drop sum
drop dfipref? dfiweight?


**************************************************************************
* Country types
**************************************************************************
/* Initial country types: */
* Half of the countries start as type 1, 1/3 as type 2, 1/6 as type 3
gen type = 1
replace type = 2 if country>`nC'*0.5
replace type = 3 if country>`nC'*(5/6)


/* Evolution of country types over time: */
* generate a 0-1 random number to base transitions on.
generate double randnum = uniform()
if "`transition'" == "default" {
	/* We want to implement the following transition matrix: .85 & .10 & .05\\ .05 & .85 & .10\\ .05 & .05 & .9
	Run transition matrix to change country types over time: */
	sort country year
	forvalues j = 2/`t' {
		// type 1 becoming type 2:
		replace type=2 if randnum>0.85 & randnum<=0.95 & L.type==1 & year==`j'
		// type 1 becoming type 3:
		replace type=3 if randnum>0.95 & L.type==1 & year==`j'
		// type 1 staying type 1:
		replace type=1 if randnum<= 0.85 & L.type==1 & year==`j'
		
		// type 2 becoming type 1:
		replace type=1 if randnum<=0.05 & L.type==2 & year==`j'
		// type 2 becoming type 3
		replace type=3 if randnum>0.9 & L.type==2 & year==`j'
		// type 2 staying type 2:
		replace type=2 if randnum>0.05 & randnum<=0.9 & L.type==2 & year==`j'
		
		// type 3 becoming type 1:
		replace type=1 if randnum<=0.05 & L.type==3 & year==`j'
		// type 3 becoming type 2:
		replace type=2 if randnum>0.05 & randnum<=0.1 & L.type==3 & year==`j'
		// type 3 staying type 3:
		replace type=3 if randnum>0.1 & L.type==3 & year==`j'
	}
}
else if "`transition'" == "fixedtypes" {
	// do nothing; this keeps types fixed over time
}
else if "`transition'" == "oneperiodchange" {
	/* Use same probabilities as in default set-up but now these are the probabilities in each period that you are of a certain type.
	In other words, deviations from core type from period 2 onwards are temporary; keep reverting back to core type. */
	sort country year
	gen coretype = type
	// type 1 becoming type 2:
	replace type=2 if randnum>0.85 & randnum<=0.95 & coretype==1 & year > 1
	// type 1 becoming type 3:
	replace type=3 if randnum>0.95 & coretype==1 & year > 1
	
	// type 2 becoming type 1:
	replace type=1 if randnum<=0.05 & coretype==2 & year > 1
	// type 2 becoming type 3
	replace type=3 if randnum>0.9 & coretype==2 & year > 1
	
	// type 3 becoming type 1:
	replace type=1 if randnum<=0.05 & coretype==3 & year > 1
	// type 3 becoming type 2:
	replace type=2 if randnum>0.05 & randnum<=0.1 & coretype==3 & year > 1
	
	drop coretype
}
else if "`transition'" == "fewertransitions" {
	/* We now want to implement the following transition matrix: .97 & .02 & .01\\ .01 & .97 & .02\\ .01 & .01 & .98
	Run transition matrix to change country types over time: */
	sort country year
	forvalues j = 2/`t' {
		// type 1 becoming type 2:
		replace type=2 if randnum>0.97 & randnum<=0.99 & L.type==1 & year==`j'
		// type 1 becoming type 3:
		replace type=3 if randnum>0.99 & L.type==1 & year==`j'
		// type 1 staying type 1:
		replace type=1 if randnum<= 0.97 & L.type==1 & year==`j'
		
		// type 2 becoming type 1:
		replace type=1 if randnum<=0.01 & L.type==2 & year==`j'
		// type 2 becoming type 3
		replace type=3 if randnum>0.98 & L.type==2 & year==`j'
		// type 2 staying type 2:
		replace type=2 if randnum>0.01 & randnum<=0.98 & L.type==2 & year==`j'
		
		// type 3 becoming type 1:
		replace type=1 if randnum<=0.01 & L.type==3 & year==`j'
		// type 3 becoming type 2:
		replace type=2 if randnum>0.01 & randnum<=0.02 & L.type==3 & year==`j'
		// type 3 staying type 3:
		replace type=3 if randnum>0.02 & L.type==3 & year==`j'
	}
}
else {
	di as error "Incorrect transition of country types specified"
}
capture drop randnum


**************************************************************************
* Budgets
**************************************************************************
/* Time-varying budget for DFI sector as a whole. */
if "`dfibudget'" == "default" {
	sort country year
	gen int db = `db' in 1/`t' // this just does the first country
	* set budget
	replace db = round(L.db*exp(`drift' + rnormal(0,`sigma'))) if year > 1 & db != .
	gegen double db2 = max(db), by(year) // this fills in the same db for each year for the other countries
	drop db
	rename db2 db
}
else if "`dfibudget'" == "capped" {
	sort country year
	gen DFIbudgetshare = runiform(0,0.1) in 1/`t'
	gen int db = DFIbudgetshare*(`nI'*`nC') in 1/`t' // this just does the first country
	* set budget
	gegen double db2 = max(db), by(year) // this fills in the same db for each year for the other countries
	drop db
	rename db2 db
}
else {
	di as error "Incorrect DFI budget specification"
}


/* Private budget. */
* Does not need to be stochastic, all that matters is whether pb is larger than bankable set of investments (as specified at the start of the code)
gen double pb = `pb'


**************************************************************************
* Individual projects
**************************************************************************
/* Data for individual projects */
expand `nI'
sort country year
* create an index for individual projects
gen i = _n
order country year i


/* Project characteristics */
gen double pc = rnormal(`mu1',`musd1')
replace pc = rnormal(`mu2',`musd2') if type==2
replace pc = rnormal(`mu3',`musd3') if type==3
label var pc "Project characteristics"


/* Expected returns */
gen double er = pc + rnormal(0,`ersd')
label var er "Expected returns"


********************************************************************
* DFI investment selection mechanism
********************************************************************
if "`mechanism'" == "random" {
	/* DFI set of eligible projects is based on expected return thresholds. 
	Within this global set, in each year the DFI randomly picks projects until its budget runs out. */


	/* Selection of investments by DFI sector as a whole. */
	* Set of investments DFIs would like to invest in:
	gen dfiset = inrange(er,`dfilo',`dfihi')

	
	/* In each year, DFI sector picks investments from this global set of investment opportunities (so not within each county) until the budget is exhausted */
	/* Put investable projects in random order and pick from the top until the budget is exhausted, 
	so that DFI sector randomly picks projects from its global set of investment opportunities until it exhausts its budget. */
	gen double randnum = uniform()
	sort year randnum // randomly sort projects within each year
	* cumulative count of number of projects in each year DFI is interested in:
	by year: gen cum_dfi = sum(dfiset) 
	* DFI sector invests if project is in eligible set and still within the budget:
	generate byte dfi_inv = dfiset == 1 & cum_dfi <= db 
	

	* Divide up these investments between multiple DFIs
	forvalues j = 1/`nD' {
		gen double dfi_inv_`j' = dfishare`j'*dfi_inv
	}


	/* Selection of investment by private investors: */
	gen psset = 0
	* DFI sector has first dibs so private sector can only pick from investments not already taken by DFI:
	replace psset = 1 if er >= `psmin' & dfi_inv == 0
	* same trick to have private investors pick randomly from global set of projects they are interested in in each year
	sort year randnum
	by year: gen cum_ps = sum(psset)
	* private sector invests if the project is in the eligible set and still within the budget:
	generate byte ps_inv = psset == 1 & cum_ps <= pb
}
else if "`mechanism'" == "pclow" {
	/* DFI set of eligible projects is based on expected return thresholds.
	Within this global set, in each year the DFI now first picks the projects with the worst project characteristics. */


	/* Selection of investments by DFI sector as a whole. */
	* Set of investments DFIs would like to invest in:
	gen dfiset = inrange(er,`dfilo',`dfihi')
	/* In each year, DFI sector picks investments from this global set of investment opportunities (so not within each country) until the budget is exhausted */
	/* Rank investment opportunities within each year by pc and pick lowest pc projects first, until budget is exhausted. */
	sort year pc
	* cumulative count of number of projects in each year DFI is interested in:
	by year: gen cum_dfi = sum(dfiset)
	* DFI invests if project is in eligible set and still within the budget:
	generate byte dfi_inv = dfiset == 1 & cum_dfi <= db


	* Divide up these investments between multiple DFIs
	forvalues j = 1/`nD' {
		gen double dfi_inv_`j' = dfishare`j'*dfi_inv
	}

	
	/* Selection of investment by private investors: */
	gen psset=0
	* DFI has first dibs so private sector can only pick from investments not already taken by DFI:
	replace psset = 1 if er >= `psmin' & dfi_inv == 0
	* have private investors pick randomly from global set of projects they are interested in in each year
	gen double randnum = uniform()
	sort year randnum
	by year: gen cum_ps = sum(psset)
	* private sector invests if the project is in the eligible set and still within the budget:
	generate byte ps_inv = psset == 1 & cum_ps <= pb 
}
else {
	di as error "Incorrect investment selection mechanism specified"
}
drop dfishare? cum_dfi psset randnum cum_ps


/*
***************************************************************
* Check percentages of bankable projects in single large draw 
***************************************************************
/* Check what percentage of investment opportunities are bankable for the private and DFI sectors. 
This is done in a single run with large nI (50000) and with thresholds for the DFI set between 0 and 2; after that this code is commented out. 
(One way to do this single run is to comment out the definition of the elusive program above and to define a set of locals instead for the parameters used in the simulation.) */
// percentage of projects bankable for private sector and DFIs, and median project characteristics by type
foreach num of numlist 1/3 {
	count if type == `num'
	scalar type`num' = r(N)
	count if er > `psmin' & type == `num'
	scalar type`num'bankpriv = r(N)
	count if dfiset == 1 & type == `num'
	scalar type`num'bankDFI = r(N)
	di "Percentage of bankable projects for type `num' for private sector is " (type`num'bankpriv/type`num')*100
	di "Percentage of bankable projects for type `num' for DFIs is " (type`num'bankDFI/type`num')*100
	su pc if type == `num' & dfiset == 1, detail
}
*/
drop dfiset


******************************************************************************
* Number of projects with expected return greater than 2 in each country-year
******************************************************************************
gen erover2 = 1 if er >= 2
bysort country year (i): gegen numprojerover2 = count(erover2)
drop erover2


******************************************************************************
* Move to dyadic format
******************************************************************************
/* create different DFIs */
sort country year
expand `nD'
sort i
gegen dfi_id = seq(), f(1) t(`nD')
label var dfi_id "DFI name"
gen double dfiinvestment = 0
label var dfiinvestment "Investments by DFIs"
forvalues j = 1/`nD' {
	replace dfiinvestment = dfi_inv_`j' if dfi_id == `j'
}
drop dfi_inv_?
order country year dfi_id i


********************************************************************************
* Collapse data to country-year-DFI format by summing over projects
********************************************************************************
sort country year dfi_id i
gcollapse (mean) type db pc numprojerover2 (sum) dfiinvestment ps_inv, by(country year dfi_id) 
label var db "DFI sector budget"
label var pc "Mean project characterstics"
/* WARNING: now ps_inv is there multiple times (once for each DFI)
so do not sum over DFIs within country-period to calculate total private investment for each country-year. */


************************************************************
* Generate supply-push instrument
************************************************************
sort country dfi_id year
// calculate budget of each DFI in each year. Summed across DFIs, this equals the variable db
bysort dfi_id year (country): gegen double DFIbudgetannual = sum(dfiinvestment)
// calculate each DFI's investments in each country as a share of each DFI's annual budget:
gen double shareRDFIannual = (dfiinvestment/DFIbudgetannual) 
/*
// double check that shares sum to 1 when summed over all recipient countries within year and DFI:
bysort dfi_id year (country): egen double sharecheck = sum(shareRDFIannual)
su sharecheck
drop sharecheck
*/
// calculate initial share by averaging annual shares over the first 5 years:
bysort country dfi_id (year): gegen double shareRDFIinittemp = mean(shareRDFIannual) if year <= 5
bysort country dfi_id (year): gegen double shareRDFIinit = max(shareRDFIinittemp)
drop shareRDFIinittemp
// generate instrument as initial share*total annual budget of DFI (sum of instrument taken over DFIs later, when collapsing)
gen double DFIinvIV = shareRDFIinit*DFIbudgetannual
*browse country year dfi_id dfiinvestment DFIbudgetannual db shareRDFIinit DFIinvIV


/* Leave one out version of the supply-push instrument */
// calculate budget of each DFI in each year minus country's own DFI investment. 
// DFIbudgetannuall1out is total amount that DFI invests in that year in other countries.
gen double DFIbudgetannuall1out = DFIbudgetannual - dfiinvestment
// generate instrument as initial share*annual budget of DFI
gen double DFIinvIVl1out = shareRDFIinit*DFIbudgetannuall1out
*browse country year dfi_id dfiinvestment DFIbudgetannual db shareRDFIinit DFIinvIV DFIbudgetannuall1out DFIinvIVl1out


********************************************************************************
* Collapse data to country-year format
********************************************************************************
sort country year dfi_id
d, fulln
*browse
gcollapse (mean) type pc ps_inv numprojerover2 db (sum) dfiinvestment DFIinvIV DFIinvIVl1out, by(country year)


********************************************************************************
* REGRESSIONS
********************************************************************************
/* Calculate total investment for regressions: */
sort country year
// to test billions to trillions strategy (conditional on starting with full additionality)
replace ps_inv = ps_inv + (`btot'-1)*dfiinvestment
// e.g. if btot = 3, then 2*DFI investment is added to private investment, so true coefficient on DFI investment becomes 3 when set-up starts from full additionality
// if btot = 1 (the default) then nothing changes.
gen double totalinv = dfiinvestment + ps_inv


/* Generate time dummies */
su year
local yearmax = r(max)
foreach num of numlist 1(1)`yearmax' {
	gen ydummy`num'=(year==`num')
}


/* add noise to pc */
replace pc = pc + rnormal(0,`noise')


if "`estimator'" == "OLSFE" {

	/* OLS */
	reg totalinv dfiinvestment ydummy*, cluster(country)
	return scalar olsbetadfi = _b[dfiinvestment]
	return scalar olssedfi = _se[dfiinvestment]
	return scalar olsr2 = e(r2)
	test dfiinvestment
	return scalar olspzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) // one-sided test of null that beta <= 0 against beta > 0 (null is zero additionality)
	return scalar olspfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) // one-sided test of null that beta >= 1 against beta < 1 (null is full additionality)
	/* OLS controlling for pc  */
	reg totalinv dfiinvestment pc ydummy*, cluster(country)
	return scalar olspcbetadfi = _b[dfiinvestment]
	return scalar olspcsedfi = _se[dfiinvestment]
	return scalar olspcr2 = e(r2)
	test dfiinvestment
	return scalar olspcpzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar olspcpfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	/* OLS with pc only (to pick up R2, with and without time dummies) */
	reg totalinv pc ydummy*
	return scalar olspcnodfir2 = e(r2)
	reg totalinv pc
	return scalar olspcnodfinotdummiesr2 = e(r2)

	/* FE */
	xtreg totalinv dfiinvestment ydummy*, fe cluster(country) 
	return scalar febetadfi = _b[dfiinvestment]
	return scalar fesedfi = _se[dfiinvestment]
	return scalar fewr2 = e(r2_w) 
	test dfiinvestment
	return scalar fepzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar fepfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	/* FE controlling for pc */
	xtreg totalinv dfiinvestment pc ydummy*, fe cluster(country) 
	return scalar fepcbetadfi = _b[dfiinvestment]
	return scalar fepcsedfi = _se[dfiinvestment]
	return scalar fepcwr2 = e(r2_w) 
	test dfiinvestment
	return scalar fepcpzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar fepcpfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	/* FE with pc only (to pick up within R2, with and without time dummies) */
	xtreg totalinv pc ydummy*, fe
	return scalar fepcnodfiwr2 = e(r2_w) 
	xtreg totalinv pc, fe
	return scalar fepcnodfinotdummieswr2 = e(r2_w) 

	/* Check omitted variable bias in OLS manually. 
	Can think of our set-up as having as the omitted variable the number of projects with expected return over 2 (which has a coefficient of 1). 
	Controlling for this variable should give us the true model (and a correct estimate of additionality). */
	/* OLS, true model without pc. */
	reg totalinv dfiinvestment ydummy* numprojerover2
	return scalar olsbetadfitrue = _b[dfiinvestment]
	return scalar olsbetanumprojerover2 = _b[numprojerover2]
	// auxiliary regression needed to calculate ommited variable bias:
	reg numprojerover2 dfiinvestment ydummy*
	return scalar olsdelta = _b[dfiinvestment]
	// components of delta
	// first partial out controls
	reg numprojerover2 ydummy*
	predict numprojerover2res, res
	reg dfiinvestment ydummy*
	predict dfiinvestmentres, res
	regress numprojerover2res dfiinvestmentres, nocons
	return scalar olsdeltapartialledout = _b[dfiinvestmentres]
	gen double dfiinvestmentressq = dfiinvestmentres^2
	su dfiinvestmentressq
	return scalar olsdeltadenom = r(sum)
	gen double cov = numprojerover2res*dfiinvestmentres
	su cov
	return scalar olsdeltanum = r(sum)
	// raw variance and covariance
	correlate numprojerover2 dfiinvestment, cov
	return scalar covnumprojver2DFIinv = r(cov_12) 
	su dfiinvestment
	return scalar varDFIinv = r(Var) 
	
	/* OLS, true model, now also controlling for pc. */
	reg totalinv dfiinvestment pc ydummy* numprojerover2
	return scalar olspcbetadfitrue = _b[dfiinvestment]
	return scalar olspcbetanumprojerover2 = _b[numprojerover2]
	return scalar olspcbetapctrue = _b[pc]
	// auxiliary regression:
	reg numprojerover2 pc dfiinvestment ydummy*
	return scalar olspcdelta = _b[dfiinvestment]
	// components of delta
	// first partial out controls
	reg numprojerover2 ydummy* pc
	predict numprojerover2respc, res
	reg dfiinvestment ydummy* pc
	predict dfiinvestmentrespc, res
	regress numprojerover2respc dfiinvestmentrespc, nocons
	return scalar olspcdeltapartialledout = _b[dfiinvestmentrespc]
	gen double dfiinvestmentrespcsq = dfiinvestmentrespc^2
	su dfiinvestmentrespcsq
	return scalar olspcdeltadenom = r(sum)
	gen double pccov = numprojerover2respc*dfiinvestmentrespc
	su pccov
	return scalar olspcdeltanum = r(sum)

	
	/* Pick up budget over time and number of projects with expected return over 2. */
	gcollapse (mean) db (sum) numprojerover2 dfiinvestment ps_inv, by(year)
	gen double dbpercnumprojerover2 = db/numprojerover2
	su dbpercnumprojerover2
	return scalar maxdbpercnumprojerover2 = r(max)
	//
	gen double invDFIvsprivate = dfiinvestment/ps_inv
	su invDFIvsprivate, detail
	return scalar medinvDFIvsprivate = r(p50)
	
	
}
else if "`estimator'" == "GMM" {

	gen double ltotalinv = L.totalinv
	
	/* We use a counter below for the different models:
	- #1: system GMM with dfiinvestment (and pc) endogenous, static model; 1 lagged level used as instruments in the differenced equations.
	- #2: same as #1 but with LDV.
	- #3: same as #1 but with 4 lags as instruments. 
	- #4: same as #1 but with LDV and with 4 lags as instrument. 
	- #5: building on #1, difference GMM.
	- #6: building on #1, IV on differenced equation.
	- #7: building on #1, GMM on levels equation.
	- #8: building on #1, IV on levels equation. */
	
	
	/* 1: SYSTEM GMM, STATIC MODEL, 1 LAG AS INSTRUMENT. */
	local counter 1
	xtabond2 totalinv dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small
	return scalar betadfi`counter' = _b[dfiinvestment]
	return scalar numinstr`counter' = e(j)
	return scalar m1p`counter' = e(ar1p) 
	return scalar m1`counter' = e(ar1)
	return scalar m2p`counter' = e(ar2p) 
	return scalar hansp`counter' = e(hansenp) 
	return scalar sargp`counter' = e(sarganp) 
	matrix define A`counter' = e(diffsargan)
	matrix list A`counter'
	test dfiinvestment
	return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
	return scalar condFDFI`counter' = e(hansen)
	return scalar condFDFIp`counter' = e(hansenp) 
	
	// now with pc
	xtabond2 totalinv dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse split) gmm(pc, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small
	return scalar betadfiwpc`counter' = _b[dfiinvestment]
	return scalar numinstrwpc`counter' = e(j)
	return scalar m1pwpc`counter' = e(ar1p)
	return scalar m1wpc`counter' = e(ar1) 
	return scalar m2pwpc`counter' = e(ar2p) 
	return scalar hanspwpc`counter' = e(hansenp) 
	return scalar sargpwpc`counter' = e(sarganp) 
	matrix define Awpc`counter' = e(diffsargan)
	matrix list Awpc`counter'
	return scalar diffHp_DFIinvdiffwpc`counter' = Awpc`counter'[5,2] // moment conditions related to DFI investment in the differenced equations
	return scalar diffHp_DFIinvlevwpc`counter' = Awpc`counter'[5,3] // moment conditions related to DFI investment in the levels equations
	test dfiinvestment
	return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse split) gmm(pc, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
	return scalar condFDFIwpc`counter' = e(hansen)
	return scalar condFDFIpwpc`counter' = e(hansenp) 

	
	if "`extraGMM'" == "yes" {
		/* 2: SAME AS #1 BUT WITH LDV. */
		local counter 2
		xtabond2 totalinv ltotalinv dfiinvestment ydummy3-ydummy20, gmm(totalinv, laglimits(2 2) collapse) gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
		iv(ydummy3-ydummy20, eq(level)) robust small
		return scalar betadfi`counter' = _b[dfiinvestment]
		return scalar betalagdep`counter' = _b[ltotalinv]
		return scalar numinstr`counter' = e(j)
		return scalar m1p`counter' = e(ar1p) 
		return scalar m1`counter' = e(ar1) 
		return scalar m2p`counter' = e(ar2p) 
		return scalar hansp`counter' = e(hansenp) 
		return scalar sargp`counter' = e(sarganp) 
		matrix define A`counter' = e(diffsargan)
		matrix list A`counter'
		return scalar diffHp_DFIinvdiff`counter' = A`counter'[5,3] // moment conditions related to DFI investment in the differenced equations
		return scalar diffHp_DFIinvlev`counter' = A`counter'[5,4] // moment conditions related to DFI investment in the levels equations
		test dfiinvestment
		return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment ltotalinv ydummy3-ydummy20, gmm(totalinv, laglimits(2 2) collapse) gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
		iv(ydummy3-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFI`counter' = e(hansen) 
		return scalar condFDFIp`counter' = e(hansenp) 
		
		// now with pc
		xtabond2 totalinv ltotalinv dfiinvestment pc ydummy3-ydummy20, gmm(totalinv, laglimits(2 2) collapse) gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
		gmm(pc, laglimits(2 2) collapse) iv(ydummy3-ydummy20, eq(level)) robust small
		return scalar betadfiwpc`counter' = _b[dfiinvestment]
		return scalar betalagdepwpc`counter' = _b[ltotalinv]
		return scalar numinstrwpc`counter' = e(j)
		return scalar m1pwpc`counter' = e(ar1p) 
		return scalar m1wpc`counter' = e(ar1) 
		return scalar m2pwpc`counter' = e(ar2p) 
		return scalar hanspwpc`counter' = e(hansenp) 
		return scalar sargpwpc`counter' = e(sarganp) 
		matrix define Awpc`counter' = e(diffsargan)
		matrix list Awpc`counter'
		return scalar diffHp_lagdepwpc`counter' = Awpc`counter'[5,2] // moment conditions for lag dep var
		return scalar diffHp_DFIinvdiffwpc`counter' = Awpc`counter'[5,3] // moment conditions related to DFI investment in the differenced equations
		return scalar diffHp_DFIinvlevwpc`counter' = Awpc`counter'[5,4] // moment conditions related to DFI investment in the levels equations
		return scalar diffHp_pc`counter' = Awpc`counter'[5,5] // moment conditions related to pc
		test dfiinvestment
		return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment ltotalinv pc ydummy3-ydummy20, gmm(totalinv, laglimits(2 2) collapse) gmm(dfiinvestment, laglimits(2 2) collapse split) /// 
		gmm(pc, laglimits(2 2) collapse) iv(ydummy3-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFIwpc`counter' = e(hansen) 
		return scalar condFDFIpwpc`counter' = e(hansenp)

		
		/* 3: SAME AS #1 (STATIC) BUT WITH 4 LAGS AS INSTRUMENTS FOR THE DIFFERENCED EQUATION. */
		local counter 3
		xtabond2 totalinv dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		iv(ydummy2-ydummy20, eq(level)) robust small	
		return scalar betadfi`counter' = _b[dfiinvestment]
		return scalar numinstr`counter' = e(j)
		return scalar m1p`counter' = e(ar1p) 
		return scalar m1`counter' = e(ar1) 
		return scalar m2p`counter' = e(ar2p) 
		return scalar hansp`counter' = e(hansenp)
		return scalar sargp`counter' = e(sarganp)
		matrix define A`counter' = e(diffsargan)
		matrix list A`counter'
		return scalar diffHp_lev`counter' = A`counter'[5,1] // all moment conditions related to the equations in levels
		return scalar diffHp_DFIinvlev`counter' = A`counter'[5,3] // moment conditions related to DFI investment in the levels equations
		test dfiinvestment
		return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFI`counter' = e(hansen)
		return scalar condFDFIp`counter' = e(hansenp) 
		
		// now with pc
		xtabond2 totalinv dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 5) collapse split) gmm(pc, laglimits(2 5) collapse) /// 
		iv(ydummy2-ydummy20, eq(level)) robust small
		return scalar betadfiwpc`counter' = _b[dfiinvestment]
		return scalar numinstrwpc`counter' = e(j)
		return scalar m1pwpc`counter' = e(ar1p)
		return scalar m1wpc`counter' = e(ar1) 
		return scalar m2pwpc`counter' = e(ar2p) 
		return scalar hanspwpc`counter' = e(hansenp)
		return scalar sargpwpc`counter' = e(sarganp)
		matrix define Awpc`counter' = e(diffsargan)
		matrix list Awpc`counter'
		return scalar diffHp_levwpc`counter' = Awpc`counter'[5,1] // all moment conditions related to the equations in levels
		return scalar diffHp_DFIinvdiffwpc`counter' = Awpc`counter'[5,2] // moment conditions related to DFI investment in the differenced equations
		return scalar diffHp_DFIinvlevwpc`counter' = Awpc`counter'[5,3] // moment conditions related to DFI investment in the levels equations
		return scalar diffHp_pc`counter' = Awpc`counter'[5,4] // moment conditions related to pc
		test dfiinvestment
		return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 5) collapse split) gmm(pc, laglimits(2 5) collapse) /// 
		iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFIwpc`counter' = e(hansen)
		return scalar condFDFIpwpc`counter' = e(hansenp)


		/* 4: SYSTEM GMM WITH LDV, 4 LAGS AS INSTRUMENTS IN DIFFERENCED EQUATION. */
		local counter 4
		xtabond2 totalinv ltotalinv dfiinvestment ydummy3-ydummy20, gmm(totalinv, laglimits(2 5) collapse) gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		iv(ydummy3-ydummy20, eq(level)) robust small
		return scalar betadfi`counter' = _b[dfiinvestment]
		return scalar betalagdep`counter' = _b[ltotalinv]
		return scalar numinstr`counter' = e(j)
		return scalar m1p`counter' = e(ar1p) 
		return scalar m1`counter' = e(ar1) 
		return scalar m2p`counter' = e(ar2p) 
		return scalar hansp`counter' = e(hansenp) 
		return scalar sargp`counter' = e(sarganp) 
		matrix define A`counter' = e(diffsargan)
		matrix list A`counter'
		return scalar diffHp_lev`counter' = A`counter'[5,1] // all moment conditions related to the equations in levels
		return scalar diffHp_lagdep`counter' = A`counter'[5,2] // moment conditions for lag dep
		return scalar diffHp_DFIinvdiff`counter' = A`counter'[5,3] // moment conditions related to DFI investment in the differenced equations
		return scalar diffHp_DFIinvlev`counter' = A`counter'[5,4] // moment conditions related to DFI investment in the levels equations
		test dfiinvestment
		return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment])
		return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment ltotalinv ydummy3-ydummy20, gmm(totalinv, laglimits(2 5) collapse) gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		iv(ydummy3-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFI`counter' = e(hansen)
		return scalar condFDFIp`counter' = e(hansenp)

		// Now with pc
		xtabond2 totalinv ltotalinv dfiinvestment pc ydummy3-ydummy20, gmm(totalinv, laglimits(2 5) collapse) gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		gmm(pc, laglimits(2 5) collapse) iv(ydummy3-ydummy20, eq(level)) robust small
		return scalar betadfiwpc`counter' = _b[dfiinvestment]
		return scalar betalagdepwpc`counter' = _b[ltotalinv]
		return scalar numinstrwpc`counter' = e(j)
		return scalar m1pwpc`counter' = e(ar1p) 
		return scalar m1wpc`counter' = e(ar1) 
		return scalar m2pwpc`counter' = e(ar2p) 
		return scalar hanspwpc`counter' = e(hansenp) 
		return scalar sargpwpc`counter' = e(sarganp) 
		matrix define Awpc`counter' = e(diffsargan)
		matrix list Awpc`counter'
		return scalar diffHp_levwpc`counter' = Awpc`counter'[5,1] // all moment conditions related to the equations in levels
		return scalar diffHp_lagdepwpc`counter' = Awpc`counter'[5,2] // moment conditions for lag dep
		return scalar diffHp_DFIinvdiffwpc`counter' = Awpc`counter'[5,3] // moment conditions related to DFI investment in the differenced equations
		return scalar diffHp_DFIinvlevwpc`counter' = Awpc`counter'[5,4] // moment conditions related to DFI investment in the levels equations
		return scalar diffHp_pc`counter' = Awpc`counter'[5,5] // moment conditions related to pc
		test dfiinvestment
		return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
		// Windmeijer underidentification test for DFI investment:
		xtabond2 dfiinvestment ltotalinv pc ydummy3-ydummy20, gmm(totalinv, laglimits(2 5) collapse) gmm(dfiinvestment, laglimits(2 5) collapse split) /// 
		gmm(pc, laglimits(2 5) collapse) iv(ydummy3-ydummy20, eq(level)) robust small nodiffsargan
		return scalar condFDFIwpc`counter' = e(hansen)
		return scalar condFDFIpwpc`counter' = e(hansenp)
	}	
	else if "`extraGMM'" == "no" {
		/* MODEL 2 WHEN IT IS NOT RUN: */
		local counter 2
		return scalar betadfi`counter' = 0
		return scalar betalagdep`counter' = 0
		return scalar numinstr`counter' = 0
		return scalar m1p`counter' = 0 
		return scalar m1`counter' = 0 
		return scalar m2p`counter' = 0 
		return scalar hansp`counter' = 0
		return scalar sargp`counter' = 0
		return scalar diffHp_DFIinvdiff`counter' = 0
		return scalar diffHp_DFIinvlev`counter' = 0
		return scalar pzeroadd`counter' = 0 
		return scalar pfulladd`counter' = 0
		return scalar condFDFI`counter' = 0 
		return scalar condFDFIp`counter' = 0 
		// now with pc
		return scalar betadfiwpc`counter' = 0
		return scalar betalagdepwpc`counter' = 0
		return scalar numinstrwpc`counter' = 0
		return scalar m1pwpc`counter' = 0
		return scalar m1wpc`counter' = 0
		return scalar m2pwpc`counter' = 0
		return scalar hanspwpc`counter' = 0
		return scalar sargpwpc`counter' = 0
		return scalar diffHp_lagdepwpc`counter' = 0
		return scalar diffHp_DFIinvdiffwpc`counter' = 0
		return scalar diffHp_DFIinvlevwpc`counter' = 0
		return scalar diffHp_pc`counter' = 0
		return scalar pzeroaddwpc`counter' = 0
		return scalar pfulladdwpc`counter' = 0
		return scalar condFDFIwpc`counter' = 0 
		return scalar condFDFIpwpc`counter' = 0	
		
		
		/* MODEL 3 WHEN IT IS NOT RUN: */
		local counter 3
		return scalar betadfi`counter' = 0
		return scalar numinstr`counter' = 0
		return scalar m1p`counter' = 0
		return scalar m1`counter' = 0
		return scalar m2p`counter' = 0
		return scalar hansp`counter' = 0
		return scalar sargp`counter' = 0
		return scalar diffHp_lev`counter' = 0
		return scalar diffHp_DFIinvlev`counter' = 0
		return scalar pzeroadd`counter' = 0
		return scalar pfulladd`counter' = 0
		return scalar condFDFI`counter' = 0
		return scalar condFDFIp`counter' = 0
		// now with pc
		return scalar betadfiwpc`counter' = 0
		return scalar numinstrwpc`counter' = 0
		return scalar m1pwpc`counter' = 0
		return scalar m1wpc`counter' = 0
		return scalar m2pwpc`counter' = 0
		return scalar hanspwpc`counter' = 0
		return scalar sargpwpc`counter' = 0
		return scalar diffHp_levwpc`counter' = 0
		return scalar diffHp_DFIinvdiffwpc`counter' = 0
		return scalar diffHp_DFIinvlevwpc`counter' = 0
		return scalar diffHp_pc`counter' = 0
		return scalar pzeroaddwpc`counter' = 0
		return scalar pfulladdwpc`counter' = 0
		return scalar condFDFIwpc`counter' = 0
		return scalar condFDFIpwpc`counter' = 0
		
		
		/* MODEL 4 WHEN IT IS NOT RUN: */
		local counter 4
		return scalar betadfi`counter' = 0
		return scalar betalagdep`counter' = 0
		return scalar numinstr`counter' = 0
		return scalar m1p`counter' = 0
		return scalar m1`counter' = 0
		return scalar m2p`counter' = 0
		return scalar hansp`counter' = 0
		return scalar sargp`counter' = 0
		return scalar diffHp_lev`counter' = 0
		return scalar diffHp_lagdep`counter' = 0
		return scalar diffHp_DFIinvdiff`counter' = 0
		return scalar diffHp_DFIinvlev`counter' = 0
		return scalar pzeroadd`counter' = 0
		return scalar pfulladd`counter' = 0
		return scalar condFDFI`counter' = 0
		return scalar condFDFIp`counter' = 0
		// Now with pc
		return scalar betadfiwpc`counter' = 0
		return scalar betalagdepwpc`counter' = 0
		return scalar numinstrwpc`counter' = 0
		return scalar m1pwpc`counter' = 0
		return scalar m1wpc`counter' = 0
		return scalar m2pwpc`counter' = 0
		return scalar hanspwpc`counter' = 0
		return scalar sargpwpc`counter' = 0
		return scalar diffHp_levwpc`counter' = 0
		return scalar diffHp_lagdepwpc`counter' = 0
		return scalar diffHp_DFIinvdiffwpc`counter' = 0
		return scalar diffHp_DFIinvlevwpc`counter' = 0
		return scalar diffHp_pc`counter' = 0
		return scalar pzeroaddwpc`counter' = 0
		return scalar pfulladdwpc`counter' = 0
		return scalar condFDFIwpc`counter' = 0
		return scalar condFDFIpwpc`counter' = 0
	}	
	else {
		di as error "Option extraGMM incorrectly specified"
	}
	
	
	/* 5: BUILDING ON #1, DIFFERENCE GMM VERSION. */
	local counter 5
	xtabond2 totalinv dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20) robust small noleveleq
	return scalar betadfi`counter' = _b[dfiinvestment]
	return scalar numinstr`counter' = e(j)
	return scalar m1p`counter' = e(ar1p) 
	return scalar m1`counter' = e(ar1) 
	return scalar m2p`counter' = e(ar2p) 
	test dfiinvestment
	return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20) robust small noleveleq nodiffsargan
	return scalar condFDFI`counter' = e(hansen)
	return scalar condFDFIp`counter' = e(hansenp) 

	// now with pc
	xtabond2 totalinv dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse) gmm(pc, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20) robust small noleveleq
	return scalar betadfiwpc`counter' = _b[dfiinvestment]
	return scalar numinstrwpc`counter' = e(j)
	return scalar m1pwpc`counter' = e(ar1p)
	return scalar m1wpc`counter' = e(ar1)
	return scalar m2pwpc`counter' = e(ar2p)
	test dfiinvestment
	return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(2 2) collapse) gmm(pc, laglimits(2 2) collapse) /// 
	iv(ydummy2-ydummy20) robust small noleveleq nodiffsargan
	return scalar condFDFIwpc`counter' = e(hansen)
	return scalar condFDFIpwpc`counter' = e(hansenp) 

	
	/* 6: BUILDING ON #1, IV ON THE DIFFERENCED EQUATION (ANDERSON-HSIAO). */
	local counter 6
	ivreg2 d.totalinv (d.dfiinvestment = L2.dfiinvestment) ydummy3-ydummy20, cluster(country) first savefirst small nocons
	capture drop sample
	gen sample = e(sample)
	return scalar betadfi`counter' = _b[d.dfiinvestment]
	return scalar FSrobustF`counter' = e(widstat)
	test d.dfiinvestment
	return scalar pzeroadd`counter' = ttail(r(df_r), _b[d.dfiinvestment]/_se[d.dfiinvestment])
	return scalar pfulladd`counter' = ttail(r(df_r), -(_b[d.dfiinvestment]-1)/_se[d.dfiinvestment])
	// pick up first stage coefficient
	reg d.dfiinvestment L2.dfiinvestment ydummy3-ydummy20 if sample == 1, nocons
	return scalar FSbetaL2DFI`counter' = _b[L2.dfiinvestment]
	// same but with OV as dependent variable
	reg d.numprojerover2 L2.dfiinvestment ydummy3-ydummy20 if sample == 1, nocons
	return scalar OVbetaL2DFI`counter' = _b[L2.dfiinvestment]

	
	/* Investigate bias/violations of moment conditions. */
	// First partial out variables
	reg d.numprojerover2 ydummy3-ydummy20 if sample == 1, nocons
	capture drop dnumprojerover2res
	predict dnumprojerover2res if sample == 1, res
	reg d.dfiinvestment ydummy3-ydummy20 if sample == 1, nocons
	capture drop ddfiinvestmentres
	predict ddfiinvestmentres if sample == 1, res
	reg L2.dfiinvestment ydummy3-ydummy20 if sample == 1, nocons
	capture drop L2dfiinvestmentres
	predict L2dfiinvestmentres if sample == 1, res
	correlate L2dfiinvestmentres dnumprojerover2res
	return scalar corrL2DFIu`counter' = r(rho) 
	correlate L2dfiinvestmentres ddfiinvestmentres
	return scalar corrL2DFIx`counter' = r(rho)
	su dnumprojerover2res
	return scalar sdu`counter' = r(sd)
	su ddfiinvestmentres
	return scalar sdx`counter' = r(sd)

	// now with pc
	ivreg2 d.totalinv (d.dfiinvestment d.pc = L2.dfiinvestment L2.pc) ydummy3-ydummy20, cluster(country) first savefirst small nocons
	capture drop sample
	gen sample = e(sample)
	return scalar betadfiwpc`counter' = _b[d.dfiinvestment]
	return scalar FSrobustFwpc`counter' = e(widstat)
	test d.dfiinvestment
	return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[d.dfiinvestment]/_se[d.dfiinvestment])
	return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[d.dfiinvestment]-1)/_se[d.dfiinvestment])
	// pick up first stage coefficients
	reg d.dfiinvestment L2.dfiinvestment L2.pc ydummy3-ydummy20 if sample == 1, nocons
	return scalar FSbetaL2DFIwpc`counter' = _b[L2.dfiinvestment]
	return scalar FSbetaL2pc`counter' = _b[L2.pc]
	// same but with OV as dependent variable
	reg d.numprojerover2 L2.dfiinvestment L2.pc ydummy3-ydummy20 if sample == 1, nocons
	return scalar OVbetaL2DFIwpc`counter' = _b[L2.dfiinvestment]
	return scalar OVbetaL2pc`counter' = _b[L2.pc]
	
	
	/* 7: BUILDING ON #1, LEVEL GMM VERSION. */
	local counter 7
	xtabond2 totalinv dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(1 1) collapse eq(level)) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small
	return scalar betadfi`counter' = _b[dfiinvestment]
	return scalar numinstr`counter' = e(j)
	return scalar m1p`counter' = e(ar1p)
	return scalar m1`counter' = e(ar1)
	return scalar m2p`counter' = e(ar2p)
	test dfiinvestment
	return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment ydummy2-ydummy20, gmm(dfiinvestment, laglimits(1 1) collapse eq(level)) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
	return scalar condFDFI`counter' = e(hansen)
	return scalar condFDFIp`counter' = e(hansenp) 
	
	// now with pc
	xtabond2 totalinv dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(1 1) collapse eq(level)) gmm(pc, laglimits(1 1) collapse eq(level)) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small
	return scalar betadfiwpc`counter' = _b[dfiinvestment]
	return scalar numinstrwpc`counter' = e(j)
	return scalar m1pwpc`counter' = e(ar1p) 
	return scalar m1wpc`counter' = e(ar1)
	return scalar m2pwpc`counter' = e(ar2p) 
	test dfiinvestment
	return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// Windmeijer underidentification test for DFI investment:
	xtabond2 dfiinvestment pc ydummy2-ydummy20, gmm(dfiinvestment, laglimits(1 1) collapse eq(level)) gmm(pc, laglimits(1 1) collapse eq(level)) /// 
	iv(ydummy2-ydummy20, eq(level)) robust small nodiffsargan
	return scalar condFDFIwpc`counter' = e(hansen)
	return scalar condFDFIpwpc`counter' = e(hansenp) 


	/* 8: BUILDING ON #1, IV ON THE LEVEL EQUATION. */
	local counter 8
	ivreg2 totalinv (dfiinvestment = DL.dfiinvestment) ydummy4-ydummy20, cluster(country) first savefirst small
	capture drop sample
	gen sample = e(sample)
	return scalar betadfi`counter' = _b[dfiinvestment]
	return scalar FSrobustF`counter' = e(widstat)
	test dfiinvestment
	return scalar pzeroadd`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment])
	return scalar pfulladd`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// pick up first stage coefficient
	reg dfiinvestment DL.dfiinvestment ydummy4-ydummy20 if sample == 1
	return scalar FSbetaDLDFI`counter' = _b[DL.dfiinvestment]
	// same but with OV as dependent variable
	reg numprojerover2 DL.dfiinvestment ydummy4-ydummy20 if sample == 1
	return scalar OVbetaDLDFI`counter' = _b[DL.dfiinvestment]

	/* Investigate bias/violations of moment conditions. */
	// First partial out variables
	reg numprojerover2 ydummy4-ydummy20 if sample == 1
	capture drop numprojerover2res
	predict numprojerover2res if sample == 1, res
	reg dfiinvestment ydummy4-ydummy20 if sample == 1
	capture drop dfiinvestmentres
	predict dfiinvestmentres if sample == 1, res
	reg DL.dfiinvestment ydummy4-ydummy20 if sample == 1
	capture drop DLdfiinvestmentres
	predict DLdfiinvestmentres if sample == 1, res
	correlate DLdfiinvestmentres numprojerover2res
	return scalar corrDLDFIu`counter' = r(rho) 
	correlate DLdfiinvestmentres dfiinvestmentres
	return scalar corrDLDFIx`counter' = r(rho)
	su numprojerover2res
	return scalar sdu`counter' = r(sd)
	su dfiinvestmentres
	return scalar sdx`counter' = r(sd)

	// now with pc
	ivreg2 totalinv (dfiinvestment pc = DL.dfiinvestment DL.pc) ydummy4-ydummy20, cluster(country) first savefirst small
	capture drop sample
	gen sample = e(sample)
	return scalar betadfiwpc`counter' = _b[dfiinvestment]
	return scalar FSrobustFwpc`counter' = e(widstat)
	test dfiinvestment
	return scalar pzeroaddwpc`counter' = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment])
	return scalar pfulladdwpc`counter' = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment])
	// pick up first stage coefficient
	reg dfiinvestment DL.dfiinvestment DL.pc ydummy4-ydummy20 if sample == 1
	return scalar FSbetaDLDFIwpc`counter' = _b[DL.dfiinvestment]
	return scalar FSbetaDLpc`counter' = _b[DL.pc]
	// same but with OV as dependent variable
	reg numprojerover2 DL.dfiinvestment DL.pc ydummy4-ydummy20 if sample == 1
	return scalar OVbetaDLDFIwpc`counter' = _b[DL.dfiinvestment]
	return scalar OVbetaDLpc`counter' = _b[DL.pc]


	/* Pick up budget over time and number of projects with expected return over 2. */
	gcollapse (mean) db (sum) numprojerover2 dfiinvestment ps_inv, by(year)
	gen double dbpercnumprojerover2 = db/numprojerover2
	su dbpercnumprojerover2
	return scalar maxdbpercnumprojerover2 = r(max)
	su db if year == 20
	return scalar finaldb = r(mean)
	su numprojerover2 if year == 20
	return scalar finalnumprojerover2 = r(mean)
	//
	gen double invDFIvsprivate = dfiinvestment/ps_inv
	su invDFIvsprivate, detail
	return scalar medinvDFIvsprivate = r(p50)

}
else if "`estimator'" == "IV" {

	bysort country (year): gen starttype = type[1]
	foreach type of numlist 1/3 {
		su totalinv if starttype == `type' & year == 1, detail
		return scalar type`type'invy1 = r(p50)
		su dfiinvestment if starttype == `type' & year == 1, detail
		return scalar type`type'dfiinvy1 = r(p50)
		su DFIinvIV if starttype == `type' & year == 1, detail
		return scalar type`type'dfiinvIVy1 = r(p50)
		//
		su totalinv if starttype == `type' & year == 20, detail
		return scalar type`type'invy20 = r(p50)
		su dfiinvestment if starttype == `type' & year == 20, detail
		return scalar type`type'dfiinvy20 = r(p50)
		su DFIinvIV if starttype == `type' & year == 20, detail
		return scalar type`type'dfiinvIVy20 = r(p50)
	}
	drop starttype

	
	/* We want to run IV regressions excluding the first 5 years, so best to drop them here already. */
	drop if year <= 5
	drop ydummy1 ydummy2 ydummy3 ydummy4 ydummy5

	
	/* FE IV regressions */
	xtivreg2 totalinv (dfiinvestment = DFIinvIV) ydummy*, fe cluster(country) first savefirst small
	return scalar ivbetadfi =_b[dfiinvestment]
	return scalar ivsedfi =_se[dfiinvestment]
	return scalar underid = e(idstat)
	return scalar underidp = e(idp)
	return scalar weakid = e(widstat)
	test dfiinvestment
	return scalar ivpzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar ivpfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	// pick up first stage coefficient
	xtreg dfiinvestment DFIinvIV ydummy*, fe cluster(country) 
	return scalar ivFSbeta = _b[DFIinvIV]
	return scalar ivFSse = _se[DFIinvIV]
	// pick up reduced form coefficient
	xtreg totalinv DFIinvIV ydummy*, fe cluster(country) 
	return scalar ivRFbeta = _b[DFIinvIV]
	return scalar ivRFse = _se[DFIinvIV]
	/* FE equivalent (without instrument) */
	xtreg totalinv dfiinvestment ydummy*, fe cluster(country)
	return scalar febetadfi =_b[dfiinvestment]
	
		
	/* Same but controlling for pc */
	xtivreg2 totalinv (dfiinvestment = DFIinvIV) pc ydummy*, fe cluster(country) first savefirst small
	return scalar ivpcbetadfi =_b[dfiinvestment]
	return scalar ivpcsedfi =_se[dfiinvestment]
	return scalar underidpc = e(idstat)
	return scalar underidppc = e(idp)
	return scalar weakidpc = e(widstat)
	test dfiinvestment
	return scalar ivpcpzeroadd = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar ivpcpfulladd = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	// pick up first stage coefficient
	xtreg dfiinvestment DFIinvIV pc ydummy*, fe cluster(country)
	return scalar ivpcFSbeta = _b[DFIinvIV]
	return scalar ivpcFSse = _se[DFIinvIV]
	// pick up reduced form
	xtreg totalinv DFIinvIV pc ydummy*, fe cluster(country)
	return scalar ivpcRFbeta = _b[DFIinvIV]
	return scalar ivpcRFse = _se[DFIinvIV]
	/* FE equivalent (without instrument) */
	xtreg totalinv dfiinvestment pc ydummy*, fe cluster(country)
	return scalar fepcbetadfi =_b[dfiinvestment]
		
	
	/* Investigate bias. */
	/* FE IV, true model without pc */
	xtivreg2 totalinv (dfiinvestment = DFIinvIV) numprojerover2 ydummy*, fe cluster(country) first savefirst small
	return scalar ivbetadfitrue =_b[dfiinvestment]
	return scalar ivbetanumprojerover2 = _b[numprojerover2]
	// calculate bias:
	// first partial out variables
	xtreg numprojerover2 ydummy*, fe
	predict numprojerover2res, e
	xtreg dfiinvestment ydummy*, fe
	predict dfiinvestmentres, e
	xtreg DFIinvIV ydummy*, fe
	predict DFIinvIVres, e
	correlate DFIinvIVres numprojerover2res, cov
	return scalar covzu = r(cov_12) 
	correlate DFIinvIVres dfiinvestmentres, cov
	return scalar covzx = r(cov_12)
	return scalar varz = r(Var_1)
	correlate DFIinvIVres numprojerover2res
	return scalar corrzu = r(rho) 
	correlate DFIinvIVres dfiinvestmentres
	return scalar corrzx = r(rho)
	su numprojerover2res
	return scalar sdu = r(sd)
	su dfiinvestmentres
	return scalar sdx = r(sd)
	
		
	/* FE IV, true model with pc */
	xtivreg2 totalinv (dfiinvestment = DFIinvIV) numprojerover2 pc ydummy*, fe cluster(country) first savefirst small
	return scalar ivpcbetadfitrue =_b[dfiinvestment]
	return scalar ivpcbetanumprojerover2 = _b[numprojerover2]
	return scalar ivpcbetapctrue = _b[pc]
	// calculate bias:
	// first partial out variables
	xtreg numprojerover2 ydummy* pc, fe
	predict numprojerover2respc, e
	xtreg dfiinvestment ydummy* pc, fe
	predict dfiinvestmentrespc, e
	xtreg DFIinvIV ydummy* pc, fe
	predict DFIinvIVrespc, e
	correlate DFIinvIVrespc numprojerover2respc, cov
	return scalar pccovzu = r(cov_12)
	correlate DFIinvIVrespc dfiinvestmentrespc, cov
	return scalar pccovzx = r(cov_12)
	return scalar pcvarz = r(Var_1)
	correlate DFIinvIVrespc numprojerover2respc
	return scalar pccorrzu = r(rho) 
	correlate DFIinvIVrespc dfiinvestmentrespc
	return scalar pccorrzx = r(rho)
	su numprojerover2respc
	return scalar pcsdu = r(sd)
	su dfiinvestmentrespc
	return scalar pcsdx = r(sd)


	/* FE IV regressions with leave one out version of instrument */
	xtivreg2 totalinv (dfiinvestment = DFIinvIVl1out) ydummy*, fe cluster(country) first savefirst small
	return scalar ivbetadfil1out =_b[dfiinvestment]
	return scalar ivsedfil1out =_se[dfiinvestment]
	return scalar underidl1out = e(idstat)
	return scalar underidpl1out = e(idp)
	return scalar weakidl1out = e(widstat)
	test dfiinvestment
	return scalar ivpzeroaddl1out = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar ivpfulladdl1out = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	/* Same but controlling for pc */
	xtivreg2 totalinv (dfiinvestment = DFIinvIVl1out) pc ydummy*, fe cluster(country) first savefirst small
	return scalar ivpcbetadfil1out =_b[dfiinvestment]
	return scalar ivpcsedfil1out =_se[dfiinvestment]
	return scalar underidpcl1out = e(idstat)
	return scalar underidppcl1out = e(idp)
	return scalar weakidpcl1out = e(widstat)
	test dfiinvestment
	return scalar ivpcpzeroaddl1out = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
	return scalar ivpcpfulladdl1out = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
	
	if `drift' != 0 {
		/* FE IV regression with country-specific trends. */
		gen trend = year - 5
		xi: xtivreg2 totalinv (dfiinvestment = DFIinvIV) i.country|trend, fe cluster(country) first savefirst small
		return scalar ivbetadfict =_b[dfiinvestment]
		return scalar ivsedfict =_se[dfiinvestment]
		return scalar underidct = e(idstat)
		return scalar underidpct = e(idp)
		return scalar weakidct = e(widstat)
		test dfiinvestment
		return scalar ivpzeroaddct = ttail(r(df_r), _b[dfiinvestment]/_se[dfiinvestment]) 
		return scalar ivpfulladdct = ttail(r(df_r), -(_b[dfiinvestment]-1)/_se[dfiinvestment]) 
		// pick up first stage coefficient
		xi: xtreg dfiinvestment DFIinvIV i.country|trend, fe cluster(country) 
		return scalar ivFSbetact = _b[DFIinvIV]
		return scalar ivFSsect = _se[DFIinvIV]
		// pick up reduced form coefficient
		xi: xtreg totalinv DFIinvIV i.country|trend, fe cluster(country) 
		return scalar ivRFbetact = _b[DFIinvIV]
		return scalar ivRFsect = _se[DFIinvIV]
	}
	else {
		return scalar ivbetadfict = 0
		return scalar ivsedfict = 0
		return scalar underidct = 0
		return scalar underidpct = 0
		return scalar weakidct = 0
		return scalar ivpzeroaddct = 0
		return scalar ivpfulladdct = 0 
		return scalar ivFSbetact = 0
		return scalar ivFSsect = 0
		return scalar ivRFbetact = 0 
		return scalar ivRFsect = 0
	}

	
	/* Pick up budget over time and number of projects with expected return over 2. */
	gcollapse (mean) db (sum) numprojerover2 dfiinvestment ps_inv, by(year)
	gen double dbpercnumprojerover2 = db/numprojerover2
	su dbpercnumprojerover2
	return scalar mindbpercnumprojerover2 = r(min)
	return scalar maxdbpercnumprojerover2 = r(max)
	su db if year == 20
	return scalar finaldb = r(mean)
	su numprojerover2 if year == 20
	return scalar finalnumprojerover2 = r(mean)
	//
	gen double invDFIvsprivate = dfiinvestment/ps_inv
	su invDFIvsprivate, detail
	return scalar medinvDFIvsprivate = r(p50)

}
else {
	di as error "Incorrect estimator specified"
}

end // end program

